######### This script will remove small patches from a species realized distribution ASCII
#########

#### Load library

library('SDMTools')

#### Directory to read files from

in.dir = '/home1/99/jc152199/skinkpatchstat/'

### Directory to write images to

#out.dir = '/home1/99/jc152199/skinkpatchstat/boxplots/'
out.dir = '/home1/99/jc152199/skinkpatchstat/histograms/'

##### Identify species

spp = c('SAPCZEC','SAPTETR','SAPBASI','GNYQUEE','CARRUBR','LAMCOGG','LAMROBE')

### Loop through each species and each statistic, producing boxplots comparing stats within species between models

for (s in spp)

	{
	
	### Read in patch stat summaries for a species
	
	acps = read.csv(paste(in.dir,s,'_microCLIM_PatchStatSummary.csv',sep=''), header=T)
	bcps = read.csv(paste(in.dir,s,'_BIOCLIM_PatchStatSummary.csv',sep=''), header=T)
	
	#### Remove first row (stats for no occurrences) and first column (row names)
	
	acps = acps[-1,-1]
	bcps = bcps[-1,-1]
	
	### Create a folder for each species
	
	dir.create(paste(out.dir,s,'/',sep=''))
	
	#### Now begin looping through individual stats
	
	for (stat in names(acps)[-1])
	
		{
		
		### Select a single statistic to work with
		
		acss = acps[,which(names(acps)==stat)]
		bcss = bcps[,which(names(bcps)==stat)]
		
		### Now calculate the variance of the two datasets
	
		acvar = var(acss)
		bcvar = var(bcss)
		
		### Calculate DOF
		
		acdof = length(acss)-1
		bcdof = length(bcss)-1
		
		#### Calculate the F-Ratio (larger variance value is always the numerator)
		
		if (acvar>bcvar)
		
			{
			
			#### Calculate the Critical F-Value
			#### This means the variance of one data set must be CritF times greater than the other before we can determine the variances are significantly different
			#### Homogeneous variance is an important assumption for the t-test
	
			CritF = qf(.975,acdof,bcdof)
		
			### Calculate the F Statistic
		
			FRatio = acvar/bcvar
			
			### Calculate the p-value for an F Statistic with x degrees of freedom
		
			Fp = 2*(1-pf(FRatio,acdof,bcdof))
			
			}
			
		if(bcvar>acvar)	
			
			{
			
			#### Calculate the Critical F-Value
			#### This means the variance of one data set must be CritF times greater than the other before we can determine the variances are significantly different
			#### Homogeneous variance is an important assumption for the t-test
	
			CritF = qf(.975,bcdof,acdof)
			
			### Calculate the F Statistic
			
			FRatio = bcvar/acvar
			
			### Calculate the p-value for an F Statistic with x degrees of freedom
		
			Fp = 2*(1-pf(FRatio,bcdof,acdof))
			
			}
		
		### If variances are significantly different, complete t-test and bind data
		
		#if(FRatio<abs(CritF))
		
			#{
			
			### t-test
			
			tt = t.test(acss,bcss)
			
			#}
			
		#if(FRatio>=abs(CritF))
		
			#{
			
			#### Perform wilctest
			
			wilc.test = wilcox.test(acss,bcss)
			
			#}
		
		### Concatenate them
		
		boxdata = c(acss,bcss)
		
		### Create a list of labels
		
		labs = c(rep('microCLIM',length(acss)),rep('BIOCLIM',length(bcss)))
		
		### Open the .png device driver
		
		#png(paste(out.dir,'/',s,'/',s,'_',stat,'_BoxPlotComparison.png',sep=''), units='cm', height=10, width=10, res=500)
		
		### Create the boxplot
		
		#boxplot(boxdata~labs, main = paste(s,'-',stat,sep=''), varwidth=TRUE, boxwex=.8)
		
		### Turn off the device driver
		
		#dev.off()
		
		### Open the .png device driver
		
		png(paste(out.dir,'/',s,'/',s,'_',stat,'_HistogramComparison.png',sep=''), units='cm', height=7, width=14, res=750)
		
		### Parameterize plot space for side-by-side plotting of histograms
		
		par(mfrow=c(1,2))
		
		#### Calculate the breaks for both boxplots
		
		breaks = seq(min(boxdata),max(boxdata),((max(boxdata)-min(boxdata))/25))
		
		#### Calculate histograms
		
		hist1 = hist(acss,breaks=c(breaks),plot=FALSE)
		hist2 = hist(bcss,breaks=c(breaks),plot=FALSE)
		
		### Determine max y-value
		
		maxy = max(c(hist1$counts,hist2$counts))
		
		### Calculate ticks for axes
		
		xticks = seq(min(boxdata),max(boxdata),((max(boxdata)-min(boxdata))/4))
		yticks = seq(0,maxy,maxy/4)
		
		### Labels for ticks
		
		xtick.labs = c(round(xticks,2))
		ytick.labs = c(round(yticks,2))
		
		### Plot histograms
		
		hist(acss,breaks=c(breaks),freq=TRUE, main='accuCLIM', xlab=stat,cex.lab=.7, cex.main=.7,plot=TRUE, ylim=c(0,maxy),axes=F)
		axis(side=1,at=xticks,labels=xtick.labs,cex.axis=.7)
		axis(side=2,at=yticks,labels=ytick.labs,cex.axis=.7)
		
		hist(bcss,breaks=c(breaks),freq=TRUE, main='BIOCLIM', xlab=paste('Wilc PVal - ',round(wilc.test$p.value,2),' T-test PVal - ',round(tt$p.value,2),sep=''),cex.lab=.7, cex.main=.7,plot=TRUE, ylim=c(0,maxy), axes=F)
		axis(side=1,at=xticks,labels=xtick.labs,cex.axis=.7)
		axis(side=2,at=yticks,labels=ytick.labs,cex.axis=.7)
		
		### Turn of device driver
		
		dev.off()
		
		### Report progress
		
		cat('\n',s,' ',stat,' Completed\n',sep='')
		
		}
		
	#### Report progress

	cat('\n',s,' All Stats Completed\n',sep='')
		
	}
	
#### Close loop
